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Abstract 

In recent years, we are seeing the formulation and use of elaborate and 
complex models in ecological studies. The questions related to the efficient, 
systematic and error-proof exploration of parameter spaces are of great im- 

CN portance to better understand, estimate confidences and make use of the 

output from these models. In this work, we investigate some of the rele- 

Q vant questions related to parameter space exploration, in particular using the 

technique known as Latin Hypercube Sampling and focusing in quantitative 

^ output analysis. We present the analysis of a structured population growth 

model and contrast our findings with results from previously used techniques, 
known as sensitivity and elasticity analyses. We also assess how are the ques- 
tions related to parameter space analysis being currently addressed in the 
ecological literature. 
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1 Introduction 



There is a growing trend in the use of mathematical modelling tools in the study 
of many areas of the biological sciences. The use of models is essencial as they 
present an opportunity to address questions that are impossible or impractical to 
answer in either in purely theoretical analyses or in field or laboratory experiments, 
and to identify the most important proccesses which should then be investigated 
by experiments. One compelling example is made by the Individual Based Models 
(IBM), which represent individuals that move and interact in space, according to 
some decision-making rules. These models permit a great level of detail and realism 
to be included, as well as linking multiple levels of complexity in a system. 

On the other hand, more realistic models employ a vast selection of input pa- 
rameters, from temperature and rainfall to metabolic and encounter rates, which 
may be difficult to accurately measure. Moreover, one may be interested in esti- 
mating how much predictions for a model fitted in one place or to one species may 
be extrapolated to different places or species. While variations in some of those 
parameters will have negligible impact on the model output, other parameters may 
profoundly impact the validity of a model's predictions, and it may be impossible to 
determine a priori which are the most important parameters. In theory, this could 
be done by evaluating the model at all possible combinations of parameters, however 
this would require a prohibitive number of model runs, specially considering that 
a single run of those models may take days to complete. Our challenge then con- 
sists in providing the best estimates for the importance of the several parameters, 
requiring the least number of model runs. 

The disciplines of uncertainty and sensivity analysis have been developed in the 
context of the physical sciences and engineering, and have been greatly developed 
in the 1980 and 1990 decades HH [33l HOI HH HH [H HH SH [19l ESI HZl HH [271 EH - 
More recently, these analyses have been successfully applied to biological models, 
in order to explore the possible outcomes from the model output, estimate their 
probability distribution and the dependency of the output on different combinations 
of parameters, and to assess which parameters require more experimental effort in 
order to be more confidently estimated. This kind of parameter space exploration 
is considered a fundamental step prior to using the model in management decisions 

One approach to the parameter space exploration, which will be described here, 
is to generate samples from the parameter space, run the model with these samples, 
and analyse the qualitative or quantitative differences in the model output. Section 
[2] will present the sampling techniques, with empashis on the Latin Hypercube 
method, while section [3] will deal with quantitative analyses. We should emphasize 
that the analyses tools are not coupled with the sampling techniques used: one 
can, in principle, use the sampling techniques described and other analyses tools, 
or apply the analyses described here to a more general class of sampling methods. 
We also present an example of the sampling and analysis used in section |4} Then, 
we briefly review some relevant research papers which have used such techniques in 
the exploration of ecological models in section [5] 



3 



1.1 Parameter spaces 

In order to better pose our questions, we need first to discuss some properties of 
the parameter space (or PS for short) of our models. 

The parameters (or inputs) are quantities X\,X2, •■• ,x m which will be used to 
run the model. In our discussion, we will assume that all the Xi are real valued. 
These quantities are unknown, and one first challenge is to determine which set of 
values better fit a model to the available data, which is the subject of linear and 
nonlinear estimation. 

However, the same model may be parametrized in different ways, as discussed 
by Ross [31]. For example, in population ecology, the logistic growth equation may 
be represented with two parameters r and K as: 

dN ( N\ 

rN 1 - — (1) 



dt V K , 

However, the same equation may be written as 

d ^L = aN + PN 2 (2) 
dt 

Here, the two parameters are a = r and (3 = —r/K. While the first equation is far 
more commonly used in the biological context, both are equivalent, and using one 
or the other model is simply a matter of choice. 

There are many other ways of writing this equation, and one of special interest 
when trying to fit real data is in terms of orthogonal polynomials, such as: 

= 9 + 9 1 (x - x) + 9 2 ((x - xf - (x~^f) (3) 

dt 

Where 9q can be calculated from 9\ and 82. This more complicated equation has 
several numerical advantages over the previous, as the parameters 9\ and 62 can be 
estimated with much more accuracy, and will not be correlated (as is the case with 
a and (5, as well as r and K). However, these parameters are hard to interpret in 
biological terms. 

These different equations illustrate the existence of interpretable (Eq. [T|, defin- 
ing, or algebraic (Eq. [2]) and computing (Eq. [3j parameters. Most of the times, it 
would be preferable to estimate the values that best fit some data by using comput- 
ing parameters, and then to transform them to interpretable parameters in order 
to present the results. 

Also, it should be mentioned that the parameter space may be constrained. This 
will have an impact on some of the available sampling and analysis techniques. The 
simplest constraint is requiring some parameter to be positive or negative. Also, 
there may be combinations of values that are meaningless. For example, if we 
model a community with N individuals and S species, the number of individuals 
and species, considered on their own, may be any positive number. However, it is 
clear that the number of species may not be bigger than the number of individuals, 
which imposes the condition S < N. This condition is called a constraint, and 
limits the values that the parameter vector may assume. 
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If we consider the m-dimensiorial space consisting on all possible combination 
of values for the parameters, our parameter space will be the subset of this space 
that respects all our constraints. For example, consider that the parameters we are 
interested are two angles of a triangle. In this case, the sum of the angles must be 
less than 180 degrees, a\ + a 2 < 180o. Clearly, this parameter space is not square, 
in the sense that, if we define the ranges of the variables a\ and a 2 independently as 
(0, 180), not all combinations of parameters will be meaningful. What can be done 
in this case is to create a new parameter 01, defined as 

* = life (4) 

This new parameter varies between and 1, and all combinations of 01,02 are 
points from our parameter space. Now, care must be exercised after applying such 
transformations in order to preserve the marginal distributions from the original 
variables, as exemplified on figure [[} 




Figure 1: a. The constrained parameter space considered, with the line representing 
a\ + a 2 = 180. The symbols represent a uniform sample taken from the space, b. 
The transformed parameter space 01,02 (see eq. [1]), showing the same sampled 
points. 

Another related concept that should not be confused with the constraints is the 
correlation between variables. For example, acidic soils are likely to have a lower 
cation exchange capacity (CEC), and more alkaline soils are likely to have a larger 
CEC [39J . Thus, those variables are correlated. Correlations have a profound impact 
on some analysis, however, they are difficult to measure, and data on correlations 
are not generally available on the literature. We will return to questions related to 



correlations in parameter spaces in section 2.2 
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1.2 Applications of parameter space exploration 

Next, we turn our attention to the kind of problems we might want to address with 
the exploration of the parameter space. First, the simplest case is asking "is there a 
region of my parameter space where condition X holds?" This condition might be, 
for example, the extinction or coexistance of species, some pattern of distribution 
or abundance of species. We also might be interested in mapping where are these 
regions. In complex models, where several different regions might exist where the 
qualitative results of the models are very different, we may ask how many of these 
regions are there, as well as map the frontiers between them. 

Another class of problems arises when the model produces some quantitative 
response, and we are interessed in determining the dependency of this response to 
the input parameters. For example, when modeling the dynamics of a population, 
we might want to know how the final population varies with each of the input 
parameters. In this context of quantitative analysis, the questions are divided in 
two classes: first, how much the variation of the input parameters is translated into 
the total variation of the results, which is the topic of uncertainty analysis, and 
second, how much of the variation in the results can be ascribed to the variation of 
each individual parameter, which is the topic of sensitivity analysis [151 EH]- We will 
present the techniques and results from both uncertainty and sensitivity analysis in 
section El 

Also, the model which we wish to analyse can be any function of the input 
parameters. In particular, there are three classes of models that can be used. First, 
the model may be a complex mathematical function (for example, Y is defined by a 
differential equation). Second, the model may be a simulation model, like an IBM. 
Third, the model may be the result of fitting a statistical model. 

All these problems may be formulated in a general way, defining some response 
from the model Y as a function of the input parameter vector x: 

Y = f (x) (5) 

In the equation [5j all the quantities are vectors, indicated by the boldface. 
Here, x = [x\, X2, ■ ■ ■ ,x m ) represent the parameters to the model f, and Y = 
[2/1,2/2, •• • , 2/n] represent the some quantitative responses from the model. In some 
sections, we will discuss the response as a single value y, without loss of generality. 

Each of the input parameters Xi is associated with a probability distribution 
Di(x), which represent our degree of knowledge about the values that Xi may assume 
(see figure [2] for examples) . 

Taken together, all the distributions Di form the joint probability distribution 
of the parameters, D(x). This function takes into account not only the individual 
distribution of each parameter, but also all the correlation terms between them. 

In very simple models, it may be possible to analitically deduce the behaviour of 
the model response taken at each point of the joint distribution of parameters. In 
the general case, however, this is impossible, and a way of investigating the model 
is to choose some points from the joint distribution and analysing the model at each 
point. Section [2] will present some strategies for choosing these points. 
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Life expectancy of a species Number of parasites per host 

(theoretical) (phenomenological) 




10 20 30 40 50 60 0.0 0.2 0.4 0.6 0.8 1.0 

x, x. 



Figure 2: Four different possibilities for chosing the distributions D{. Pannel "a" 
shows the exponential distribution of life expectancy, which can be deduced from 
theory [5J. Pannel "b" shows a gamma distribution, which can be used to model 
the number of parasites per host, but has no theoretical derivation [I]. Pannel "c" 
shows data from an empirical study on the CO2 uptake from plants [29], and pannel 
"d" shows an example where no prior information can be used. 
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2 Sampling Techniques 



There are several strategies that can be used to choose the samples from the pa- 
rameter space that will be used as input to our model of interest. Here, we will 
present some of them, along with their limitations, to justify our choice for the 



Latin Hypercube Sampling, which we will describe in section 2.1 



One way of exploring the parameter space is by discretizing every distribution 
and running the model for every possible combination of values for all parameters. 
This is called full parameter space exploration, as done in [44], and altough it 
possesses many advantages, it may become very costly in terms of computer time. 
In addition, the number of possible combinations increases exponentially with the 
number of parameter dimensions considered. 

To circumvent the exponential increase in the number of samples, it is usual to 
explore the parameter space in the following fashion: holding all but one parame- 
ter constant, we analyse how the output of a model is affected by one parameter 
dimension at a time (as done in [16]). This analysis is refered to as individual pa- 
rameter perturbance. This kind of analysis is, however, limited by the fact that the 
combinations of changed parameters may give rise to complex and unexpected be- 
haviours. Often, algorithms based on individual parameter perturbance are used as 
a first step in order to discriminate between parameters that may have a substantial 
impact on the output, and parameters that are less relevant. 

Another viable option would be to chose N random samples from the entire 
space, in order to analyse both the effect of each parameter and the combined 
effect of changing any combined number of parameters. This sampling scheme is 
called random sampling, or Monte Carlo sampling, and has been applied to many 
biological models [20]. One important feature of the Monte Carlo sampling is that 
its accuracy does not depend on the number of dimensions of the problem, as shown 
in [22]. 

Stratified sampling strategies, which are a special case of Monte Carlo sampling, 
consist in strategies for chosing these random samples while, at the same, making 
sure that every subdivision (or strata) of the distribution is well represented. As 
shown on [21], the estimatives of the statistical properties (such as the mean or the 
variance) of the model output are better represented by stratified random sampling 
than by simple random sampling (see figure [3] for examples) . As we shall see in the 
next session, the Latin Hypercube sampling is a practical and easy to understand 
stratified sampling strategy. 

Another class of Monte Carlo methods that should be mentioned here is the 
Markov Chain Monte Carlo (MCMC), which is also used on similar analyses [22J. 
This methods consists in generating a sequence of points {x^} from the parameter 
space whose distribution converges to the joint probability distribution D(x), and 
in which each sample x 1 ^ is chosen based on the previous x^ -1 ^. MCMC methods 
perform better than LHS methods for estimating the distribution of the model re- 
sponses, however, they require a number of model runs which is orders of magnitude 
higher than LHS requirements. 
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Figure 3: Illustration of four sampling methods. While the full parameter space ex- 
ploration is clearly representative of the whole space, it requires a very large number 
of samples. The individual parameter perturbation chooses samples by holding one 
parameter constant and varying the other, and clearly cannot take into account 
interactions between parameters. The random sampling uses information about the 
whole parameter space with a small number of samples, but can oversample some 
regions while undersampling others. The Latin Hypercube (section 2.1) samples all 
the intervals with equal intensity. 
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2.1 Latin Hypercube: Definition and use 

In this section, we describe the Latin Hypercube Sampling, and show how it can be 
used to efficiently solve the questions posed in section [TJ We also discuss what are 
the available methods for obtaining the LHS. 

Firstly, let us define, in the context of statistical sampling, what is a Latin 
Square: 

Definition If we divide each side in a square in N intervals, and then take samples 
from the square, the resulting square will be called Latin if and only if there is 
exactly one sample in each row and each column. 

A Latin Hypercube is simply the generalization of the Latin Square to an arbi- 
trary number of dimensions m. From these definitions, it is clear that the number 
of samples is fixed as iV We will show how to estimate the optimal N in section 



2.5, but for now it is relevant to note that N does not depend on the number of 
parameters considered. 

We will now construct the Latin Hypercube. Let's fix our attention in one 
parameter dimension i of the parameter space. The first step we should take is to 
divide the range of Xj in N equally probable intervals. In order to do so, we will 



turn our attention to the probability distribution of Xj, defined on section 1.1 as D^. 
Recall that this probability distribution must be chosen in a way that represents 
our current understanding of the biology of the given system. This function might 
be estimated by an expert in the field, it might represent a data set from field or 
laboratory work, or in some cases it may be simply the broadest possible set of 
parameters, in some cases where the actual values are unknown or experiments are 
unfeasible (see fig. [2]). 

In possession of the distribution function Di, we must sample one point from 
each equally probable interval. There are two approaches used here: it is possible 
to choose a random value from within the interval (as proposed by [21]), or instead, 
we can use the midpoint from each interval [16] . As the statistical properties of the 
generated samples are very similar, we will use the second approach here. 

The integral of the distribution function is called the cummulative distribution 
function Fi(x). This function relates the values x that the parameter may assume 
with the probability p that the parameter is less than or equal to x. We will refer to 
the inverse of the cummulative distribution function, Ff , as the quantile function 
of the parameter Xi, as it associates every probability value p in the range (0, 1) to 
the value x such that P(xi < x) = p. We divide the range (0, 1) in N intervals of 
size 1/N, and use this quantile function to determine the x values as the midpoints 
of each interval. Summarizing, we take the N points, represented as x^, k e [1 , AT] , 
from the inverse cummulative distribution Ff (x) as: 

= Fr 1 — — (6) 



N 

The samples from each dimension are subsequently shuffled, to randomize the 
order in which each value will be used (see example on figure [4]). As the samples 
come from the distributions Di, and are only reordered, their (marginal) distribution 
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Probability distribution of x 



Figure 4: Sample normal probability distribution, with 5 samples collected from 
regions with same probability and shuffled. Note that the first sample correspond 
to the second interval, the second sample correspond to the fourth interval, and so 
on. 
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will remain that of Di. However, the joint distribution of the parameters is still 
not well defined. In particular, this simple shuffling may result in some of the 
parameters to be positively or negatively correlated with each others, which might 
be undesirable. Some techniques have been developed to eliminate these correlation 
terms or to impose different correlations between the variables, and will be presented 
on section 12.21 

It should be noted that, in the mathematical literature, it is usual to refer to a 
somewhat different object as a Latin Square: this would be a square whose sides 
are divided in N intervals, and is filled with N different symbols, such that for each 
row and column there is exactly one occurrence of each symbol, as represented in 
figure [5j 




Figure 5: A stained glass window at the Caius College, Cambridge, showing a full 
Latin Square. Notice how there is only one occurence of each color in each row and 
in each column. 
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2.2 Algorithms and extensions 

As described above, the LH sampling generates an uniform distribution of samples 
in each parametric dimension. However, there is no guarantee that the correla- 
tion between two or more parameters will be zero, and the classical algorithm from 
McKay [24J, described in the previous section, usually produces correlations as high 
as 0.3 between pairs of factors, which can difficult or even compromise further analy- 
ses. In this section, we will present one algorithm designed to take into account the 
correlation between the parameter variables [16], using a single-switch-optimized 
sample reordering scheme. We will present the general case of prescribing a corre- 
lation matrix, and will also present results for the trivial case of zero correlation 
terms. Other methods have been proposed to address this problem [J71 HHl [10], in- 
cluding methods that deal with higher-order correlation terms jUJ using orthogonal 
designs. These methods, however, either impose severe restrictions on the number 
of samples that must be chosen or are too computationally intensive. 

In order to obtain the samples with prescribed correlation terms, we define as 
Cij the desired m x m correlation matrix between the variables Xj and Xj, and 
denote by C*j the current correlation between Xi and Xj. 

The next step is done iteratively for each parameter dimension, starting with the 
second one. Suppose that the method has already been applied to i — 1, 2, I — 1, 
and we will apply it to % = I. The square sum of the errors in the correlations 
between x\ and the anterior parameters is given by 

i-i 

E = Y,( c ^-°h) 2 ( 7 ) 

k=l 

Afterwards, we calculate, for each pair of values sampled from the parameter 
dimension /, what would be the error in the correlation if they were switched. 
The pair that corresponds to the greater error reduction is then switched, and the 
procedure is repeated iteratively until the error is acceptably small. 



2.3 Stochastic models 

When dealing with stochastic models, like several relevant individual based models 
(IBM), the questions presented become complicated by the fact that running the 
same model with exactly the same parameters might wield largely different results, 
both quantitative and qualitatively. In this scenario, we must be able to diferentiate 
the variation in responses due to the variation of the parameters with the variation 
in response due to stochastic effects. 

We will proceed first by defining the variation due to the input parameters as the 
epistemic uncertainty. This uncertainty arises from the fact that we do now know 
what are the correct values for a given parameter in a given natural system, and is 



related to the probability distributions Di, presented in section [L2j The variation in 
the behaviour of the model which is caused by stochastic effects is called stochastic 
uncertainty, and is inherent to the model. 

It is important to note that the two uncertainty components are impossible to 
disentangle in general stochastic models. This has prevented the general analysis of 
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such models until recently. In recent years, studies have shown that the important 
parameters and their effects can be correctly identified by running such models 
repeatedly for the same input variables and then averaging the output [35], given 
that the following conditions are respected: 

• Sample sizes should be large, relative to the aleatory uncertainty. 

• The output values should be unimodal, that is, the output values for a given 
parameter choice should be clustered around a central value. 

• The correct analysis tools should be used (as will be discussed on session [3]). 



2.4 Adaptative refinement 

After generating the N Latin Hypercube samples and running the model with the 
determined parameters, it may become necessary to increase the number of samples. 
In simple Monte Carlo simulations, increasing the number of samples can be done 
easily by generating new random points and running the model with them. However, 
to maintain the desirable properties of the Latin Hypercube, care must be taken 
to (1) keep the marginal distribution of the points equal to the parameter original 
distibution and (2) keep the correlation terms equal to zero, or equal to the ascribed 
values. 

We refer to this increase in the number of sample points as a sampling refinement, 
and as it is done on demand after examining the already generated outputs, we refer 
to it as an adaptative sampling refinement. This technique can be applied while 
estimating the Symmetrized Blest Measure of Association (SBMA) which will be 



used to determine the optimum N size (section 2.5). 

Our construction will start by defining a new sequence of values 7j for every 
parameter dimension, with 

' k - 1/6 N 



7a 



N 



•** = F-> (^) (9) 

For every k £ [1, iV], and so 7$ has 2N elements. Following this, Huntington & 
Lyrintzis' algorithm (16] is applied on the set of 7; in order to minimize the correla- 
tion between them. The extended Latin Hypercube is formed by concatenating the 
original Xi values to the 7$ values, and obtaining a cube with 3N samples defined 
as Zi, of which just 2N must be calculated. 

The extended Latin Hypercube presents the following desirable properties: 

• The expected sample mean is equal to the distribution mean for each param- 
eter: (zi) = (Di) 

• The expected sample variance is equal to the distribution variance for each 
parameter: (&(zi)) = (a(D i )) 



The correlation between any two variables is bounded by the following expres- 



sion: p\. z . < pl_ x . + p 2 
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2.5 Measuring the concordance with increasing sample size 

We will now turn our attention to the problem of determining the optimal number 
of model runs we should apply in order to provide a good estimate of which are the 
relevant parameters for a given model. One way of proceeding is by sistematically 
increasing the number N of model runs and applying any of the sensitivity analysis 
techniques, which will be discussed on the following section. If the analyses indicate 
similar results for consecutive runs, we can presume that increasing the sample size 
will not yield major changes to the results. 

All of the sensitivity analyses present us with a list of the parameters that 
have most influence in the model output. By comparing the resulting lists from 
two experiments, we can decide to stop increasing N when the lists are sufficiently 
similar. Our problem then is to determine how similar are two vectors of ranks. In 
principle, we could apply any distance function to those vectors. However, consider 
that 3 analyses indicated that the order of the most influential parameters is: 

~H1 = 1 2 3 4 5 6 
H2 = 1 2 3 6 4 5 
H3 = 2 3 1 4 5 6 

By using standard distances (like Spearman's rho or Kendall's tau), we will see the 
same difference between HI and H2 and between HI and H3. On the other hand, 
in the context of determining the most influential parameters, we would be inclined 
to see HI and H2 as more similar than any of them to H3, as the first two preserve 
the ordering of the three first parameters. 

Iman and Conover proposed a correlation coefficient for this problem called Top- 
Down Correlation Coefficient [18], which is based on Savage Scores. This coefficient, 
know as TDCC, was extensively used for sensitivity analyses [32J. Another measure 
of concordance proposed more recently is the Symmetrized Blest Measure of Asso- 
ciation (SBMA) [12]. Recent research suggests that estimates for SBMA produce a 
smaller standard error than TDCC without the assumption that there is no corre- 
lation between the variables [23]. Thus, we propose using SBMA as a measure of 
concordance between analyses from different sample sizes. Defining the ranks from 
the first sample as Ri and the ranks for the second sample as Si, the estimator for 
the SBMA is: 

,^-*ji±± + ^± r ,s,L-^) (10) 

n — 1 n A — n z — 4 \ n + 1 / 

i=l x ' 

For the techniques that may output negative values, for instance negative corre- 
lations, the SBMA must be applied on the absolute values. Otherwise, the param- 
eters which present strong negative effects will be ranked very low, and will not be 
taken into account by the SBMA. 

We will apply SBMA for the PRCC and eFAST techniques (both of which are 
discussed on section [3l on the example on section |4| 
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Box 1: Parameter space sampling 

In order to adequately explore the parameter space of a model, the following 
steps should be taken: 

• The parameters of interest must be identified, in a way that all combinations 
of parameter values are meaningful. 

• The range and distribution of each parameter should be estimated. 

• The model to be analysed must be formulated as a function of all the param- 
eters. 

• A Latin Hypercube should be generated based on the parameter distributions. 

• The correlations between parameters, when present, must be taken into ac- 
count. 

• The model must be evaluated at each combination of parameter values. If the 
model is stochastic, it should be evaluated repeatedly at each combination. 



3 Quantitative output analysis 
3.1 Uncertainty analysis 

The first question we would like to answer, in the context of quantitative analysis, 
is what is the probability distribution of the response variable y given that we know 



the join probabilities of the input parameters x (see definitions in section 1.2 ), which 
is the subject of uncertainty analysis [15] . 

This can be done by fitting a density curve to the output y or an empiric cu- 
mulative distribution function (ecdf). If there is any theoretical reason to believe 
that the distribution of y should follow one given distribution, it is possible to fit 
this function to the actual output data and estimate the distribution parameters. If 
the joint distribution of the input parameters correspond to the actual probability 
of some natural system to exhibit some given set of parameter values (as opposed 
to the case where we have no biologically relevant estimates for some parameters), 
the estimate represented by the density and ecdf functions approaches the actual 
distribution that the variable y should present in nature. This functions may be 
used, for example, to provide confidence intervals on the model responses. 

However, this is only the case when the input variables are uncorrelated or when 
enough correlation terms have been taken into account. Smith ([38]) provides an 
example where ignoring the correlation terms leads to inaccuracies on the estimation 
of confidence intervals. 

The next reasonable step is to construct and interpret scatterplots relating the 
result to each input parameter. These scatterplots may aid in the visual identi- 
fication of patterns, and altough they cannot be used to prove any relationship 
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between the model response and input, they may direct the research effort to the 
correct analyses. There are extensive reviews of the use of scatterplots to identify 
the important factors and emerging patterns in sensitivity analyses [19|. 

We will present here some quantitative analyses tools, aimed at identifying in- 
creasingly complex patterns in the model responses. It should be stressed that no 
single tool will capture all the relations between the input and output. Instead, 
several tools should be applied to any particular model. 



3.2 Sensitivity analysis 

The question of "what is the effect of some combination of parameters to the model 
output" may be answered by testing the relation between the parameters and out- 
puts. There are extensive reviews about detecting these relations after generating 
samples with Latin Hypercubes [32j [19], so we will give just a brief overview. We 
will first note that the methods used must take into account the variation of all the 
parameters. For example, instead of calculating the correlation between the result 
and some parameter, partial correlation coefficients should be used, which discount 
the effect of all other parameters. 

These relations between the result and the parameters may be classified, in order 
of increasing complexity, as: 



Linear relation, which can be tested with the Pearson partial correlation co- 
efficient, 

• Monotonic relation, which can be tested with the Spearman partial correlation 
coefficient, 

• Trends in central location, for which the Kruskal-Wallis test may be applied, 

• Trends in variability, for which the FAST method and Sobol' indices may be 
used. 



Subsections 3.2. 1| to |3.2.4 will provide some mathematical background for each 



method, and section [4] will present examples of use of those tests. We should stress 
here that the application of one method is not enough to draw conclusions about the 
relations between the input and output variables, as these techniques test different 
hypotheses, and have different statistical powers. Instead, every model should be 
analysed by a combination of techniques, preferably one for each category outlined 
here. 



3.2.1 Linear relation 

Under the hypothesis of independence between the central location and dispersion 
of the model responses, the most straightforward relationship between y and X{ is 
the linear, represented by y ~ x«. This is case if, every time x« is increased, y 
increases by aproximately the same ammount. The Pearson correlation coefficient 
is the commonly used measure to test for a linear correlation: 
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Where a a is the variance of a and a a b is the covariance between a and b. The 
correlation coefficient is a measure of the predicted change in y when Xi is changed 
one unit, relative to its standard deviations, and, as such, approaches ±1 when there 
is a strong linear relation between the variables. The square of p, usually written 
as R 2 , measures the fraction of the variance in the output that can be accounted 
for by a linear effect of Xj. Is is usual to test the significance of this linear relation 
by a t-test [TT] . 

Other than examining the individual relationships between the parameters and 
the output, we can investigate the joint effect of several x iy as y ~ X\ + x 2 + ■ ■ ■ + x m . 
In this case, the multiple R 2 represent the fraction of the variance on the output 
due to linear effects of all the Xi considered. 

However, a measure of p close to zero does not mean that no relationship exists 
between y and X{ - for instance, x 1 + y = 1, x e [—1, 1] presents p = 0, so clearly 
other methods might be needed. 

The Partial Correlation Coefficient (PCC) between Xi and y is the measure of 
the linear effect of Xi on y after the linear effects of the remaining parameters have 
been discounted. In order to calculate the PCC, first we fit a linear model of Xi as 
a function of the remaining parameters: 

Xi ~ Xi + X 2 H h Xi-i + x i+1 H hi m (12) 

A corresponding model is done with y: 

y~x 1 +x 2 -\ h + a^+i H Yx m (13) 

The PCC is calculated as the correlation between the residuals of these two 
models: 

PCC(y, x^ =p((y-y), (a* - x t )) (14) 
3.2.2 Monotonic relation 

Let us refer to each value of y as yu and each value of Xi as x^. The rank trans- 
formation of y, represented by r(^) can be found by sorting the values and 
assigning rank 1 to the smallest, 2 to the second smallest, etc, and iV to the largest. 
The rank of Xik, r(xik), can be found in a similar way. 

If there exists a strictly monotonic relation between y and Xi, that is, if every time 
Xi increases, y either always increase or always decreases by any positive ammount, it 
should be clear that the ranks of y and Xi present a linear relationship: r(y) ~ r(xi). 

The correlation between r(y) and r(xj) is called the Spearman correlation coef- 



ficient rj yXi . The same analyses presented on section 3.2.1 can also be applied for 
the rank transformed data, including significance testing and multiple regression. 

If the procedure described to calculate the PCC is followed on rank transformed 
data, that is, if y and Xi are rank transformed and fitted as linear models of the 
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remaining parameters, the correlation between the residuals is called PRCC, or Par- 
tial Rank Correlation Coefficient. This measure is a robust indicator of monotonic 
interactions between y and x iy and is subject to significance testing as described in 
|32j . This measure will perform better with increasing N. 



3.2.3 Trends in central location 

Even if the relation between y and Xi is non monotonic, it may be important and 
well-defined. The case in which y ~ xf, Xi e (—1, 1) is a common example. This 
relation may be difficult to visualize, and sometimes may not be expressed analyti- 
cally. In these cases, the Kruskal-Wallis rank sum test may be used to indicate the 
presence of such relations [19J. 

In order to perform the test, the distribution of X{ must be divided into a number 
Ntest of disjoint intervals. The model response y is then grouped with respect to 
these intervals, and the Kruskal-Wallis test is used to investigate if the y values 
have aproximately the same distribution in each of those intervals. A low p-value 
for this test indicates that the mean and median of y is likely to be different for each 
interval considered, and thus that the parameter Xi have a (possibly non monotonic) 
relationship with y. 

The number of intervals N tes t is not fixed as any "magical number", and may 
have a large impact on the test results. It is then recommended that this test 
should be repeated with different values to obtain a more comprehensive picture of 
the interactions between Xj and y (fig. [6]). 
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Figure 6: Example of application of the Kruskal-Wallis test on the same data set, 
which present a strong quadratic component, by dividing the range in 2 intervals 
(right) or 3 intervals (left). The dashed lines are the divisions between the intervals, 
and the strong horizontal lines are the sample means for each interval. 
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3.2.4 Trends in variability 

Other than the central tendency of the results, their dispersal may be dependent on 
the input parameters. A classical approach which may be used to test whether the 
dispersal of the output is related to any input parameter is to divide the distribution 
of Xi into a number N tes t of disjoint intervals and group the model response y with 
respect to these intervals, as done on the Kruskal-Wallis test. In this case, the 
ANOVA F statistic can be used to test for equality between the y conditional to 
each class [T9] . 

A similar approach, which we will use, is to employ the eFAST indices [34j 133] . 
which is a variance decomposition method based on the FAST and Sobol' indices. 
While the Sobol' indices were described in 1969, in Russian, FAST was developed 
by Cukier et al. in 1973, and both are identical in all but one computation [TJ. 

These methods estimate what fraction of the output variance can be explained 
by variation in each parameter Xi, which is called the first-order sensitivity of Xi or 
main effect of X{. The method estimates as well the fraction which is explained by 
the higher-order interactions between X{ and all other parameters. The sum of all 
terms related to Xi is called the total-order sensitivity of 'xj. 

The eFAST method estimates the main effect of each parameter by chosing a 
periodic function fi(xi) for each parameter, where the frequency (pi of each function 
is distinct, and should, in theory, be incommensurable. Each of this functions is 
sampled N s times, and a Fourier analysis is applied to the model output. The 
Fourier coefficients at each frequency (pi is related to the main effect of the variable 
Xi. The total order sensitivity of x^ is then calculated as the fraction of the variation 
which is not explained by the complimentary of Xi (that is, all parameters but this 
one). 

There are two things that should be noted here about eFAST indices. The first 
one is that the eFAST calculation does not involve the LHS sampling scheme, and 
may require more model evaluations. Also, this method produces small positive 
total-order sensitivity estimates even for parameters which do not play any role on 
the model output, as many numeric approximations are involved. 



Box 2: Uncertainty and sensitivity 

The following questions should be addressed in uncertainty and sensitivy anal- 
yses (between parentheses are the section numbers which detail each step): 



What is the range and distribution of model responses (3.1)? 



The ecdf plot of model responses suggests any theoretical distribution (3.1)? 

The scatterplots between each variable and the model response suggest any 
relationship (3.1)? 



Is there any linear or monotonic relation between the model parameters (3.2.1 
and|3T2T2l? 



Is there any trend in central location or dispersion of the response (3.2.3 and 
3~2ll)? 
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The output from the analyses is coherent between different sample sizes (2.5)? 



Case study: a density-dependent population 
matrix model 



4.1 Model description 

In this section, we demonstrate the uses and advantages of the methods outlined 
in the previous sections by performing sensitivity analyses on a density-dependent 
model of the tropical palm Euterpe edulis (commonly known as palmito jugara). All 
the data used here was extracted from Silva M atos et al. paper [37] , which compared 
a density independent matrix model of population growth with a density dependent 
model in which the recruitment of seedlings was affected by the number of seedings 
and adult trees. Silva Matos provided results, sensitivities and elasticities for the 
density independent model that can be compared to our findings, and results for 
mean and maximum values of the density dependent model - but unfortunately, 
their methods did not allow for a full sensitivity analysis of the density dependent 
model. 

We have used the R language to perform the sampling and analysis, with most 
of the code based on the "sensitivity" package, which implements PRCC analysis 
(see section 3.2.2) and eFAST analysis (see section 3.2.4), among others. We have 
also used our implementation of Huntington & Lyrintzis' algorithm to generate zero 
correlation between LHS samples (see section 2.1). All code used is freely available 
on the web. 

The models analysed are based on a Lefkovitch matrix with seven size classes. 
The matrix used on the density-independent model is 
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(15) 



Here, Pj is the probability of a tree surviving and remaining in the same class 
(statis), Gi is the probability of a tree surviving and growing to the next class, and 
Fi is the number of offspring produced per reproductive palm. 

The dominant eigenvalue of this matrix is related to the predicted population 
growth rate. We considered the dominant eigenvalue for this matrix as the model 
output, as usually done on this modelling approach. 

The density dependent model used the same matrix, but now the growth term 
of the first size class represented a decreasing function of the population density: 
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1 + aN x v P 

Here, Ni and N 7 represent the number of seedlings and adults per patch. The 
parameters G m and a represent the maximum transition rate at low densities and 
the strength of reduction in G\ with increasing seedling densities. The remaining 
parameters z and p represent the crown area of an adult tree and the plot size 
(which is fixed as 25m 2 ), and their ratio is related to the reduction of recruitment 
due to the presence of adults, due to the fact that few seedlings are able to grow 
underneath the canopy of an adult. 

As this model does not produce any static matrix, it is not meaningful to cal- 
culate any eigenvalue. Instead, the total population corresponding to the stable 
population distribution was used as model output. 

A naive approach to estimating the parameter sensitivities of this model would 
use the stasis, growth and fecundities given. However, this would yield erroneous 
results, as the probabilities of stasis and growth for a given class are not independent, 
as Pi + Gi < 1 for all classes. As discussed on section we need to use an 
alternative parametrization for this model. 

We will represent by Sj the probability of survival for each class, calculated as 
Si = Pi + Gi, and by lowercase the probability of growth, calculated as = 
(sj — Pi) I Si. Using the notation for complementary probabilities gl = 1 — g iy we can 
now write the Lefkovitch matrix as: 



A 



si ■ 9i 
si ■ 9i 















S2 ■ 9~2 
S2 ■ 92 











S3 ■ 93 
S3 ■ 93 











S4 ■ gl 
sa ■ gi 











S5 ■ 95 
S5 ■ 95 













S6 

sa 



96 
96 



F 7 






S7 



;i7) 



The models have, respectively, 14 and 16 parameters. All analyses have been 
done with mean and standard deviation calculated from Silva Matos paper, assum- 
ing a normal distribution of parameters truncated at the [0, 1] interval for proba- 
bilites, and on [0, +oo) for the other parameters. In the case of the density depen- 
dence parameters Gm, z and a, only the mean estimate was given on the paper, so 
conservative values were used for the standard deviations. 



4.2 Results 

First, we have generated Latin Hypercubes consisting of all relevant variables for 
each model. Then, the models were run for each combination of parameters. By 
using the SBMA measure of concordance (section 2.5), we have determined that 
the sample size required for the density independent model is approximately 100, 
and between 400 and 500 for the density dependent (Table [T]). 

If we presume that the data collected is representative of our knowledge about 
each of these parameters, the probability distribution of the model responses can 
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Size 


Independent 


Dependent 


1 


50-100 


0.79 


0.59 


2 


100-200 


0.90 


0.47 


3 


200-300 


0.95 


0.80 


4 


300-400 


0.96 


0.84 


5 


400-500 


0.96 


0.84 



Table 1: Comparison of PRCC analyses by sample size for both models 



be seen as the probability that the real population of palms exhibit each value 
of the model output. Figure [7] shows these distributions, which suggest that the 
population is viable for the vast majority of parameters values in the parameter 
space considered. Also, the A calculated from the density-independent model (mean 
1.22, standard deviation 0.06), is very close to the value found by Silva Matos 
(mean 1.24 ±0.06 se). Considering the density-dependent model, the median stable 
population predicted (4902 trees in each 25m 2 plot), is comparable to, altough 
higher than, the population actually measured by the study (1960 ± 560 trees per 
plot, mean and sd calculated over three years). 

We have generated scatterplots between the result from the models and each in- 
dependent parameter, in order to visually identify the relations between the inputs 



and outputs (figs. Mto 11). It is clear from these scatterplots that the fecundity 



plays a major role on the population dynamics, and may be involved in non-linear 
interactions. Also, growth probabilities (#j) have a greater impact on the model out- 
put than survival (sj) on the density-independent model. This is to be contrasted 
with Silva Matos results, which show all of the elasticities to be aproximately equal 
for all parameters. In the density- dependent model, the patterns are much more 
complex. Survival parameters seem to be more influent than growth, and the pa- 
rameters reducing the recruitment (a and z) show a clear negative effect on the 
population size. However, there is evidence now for non-linear effects of the param- 
eters, in particular si and F 7 . 

These scatterplots show very high dispersion of values, mostly due to the fact 
that all parameters are being varied between runs. In order to investigate the effect 
of each parameter on the outputs discounting the effects of the others, we analyse 



the Partial Rank Correlation Coefficient (PRCC, fig. 12). The PRCC analysis 
for the density independent model is in agreement with our previous expectations, 
with Fj being the most influential parameter, followed by growth probabilities. 
Survival probabilities follow with low correlations. The density dependent model 
presents us with some surprises, as now the survival parameter for the smallest and 
largest size classes jumped to occupy the second and third largest positive PRCCs. 
The remaining parameters follow the new parameters a and z, which are strongly 
negatively correlated with the output. 

The last analysis we present here is the Extended Fourier Amplitude Sensitivity 



Test (eFAST, fig. 13). This analysis provides an estimation of the fraction of 



variation of model output that can be explained by the individual variation of each 
parameter (which we call first-order sensitivity, or main effect), along with the total 
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Figure 7: Empirical cumulative distribution functions (ecdf) for the density inde- 
pendent (a) and density dependent (b) models of population growth. In (a), the x 
axis represents the dominant eigenvalue, and the population is viable if x > 1. In 
(b), the x axis represents the total equilibrium population, and the population is 
viable if x > 0. 
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Figure 8: Scatterplots relating the value of the input parameters of growth to the 
A calculated the output for the density independent model. Each plot shows the 
Pearson correlation coefficient and significance level: '.' p < 0.1 '*' p < 0.05 '**' p 
< 0.01 '***' p < 0.001 
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Figure 9: Scatterplots relating the value of the input parameters of survival and 
fecundity to the A calculated the output for the density independent model. Each 
plot shows the Pearson correlation coefficient and significance level: '.' p < 0.1 '*' 
p < 0.05 '**' p < 0.01 '***' p < 0.001 
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Figure 10: Scatterplots relating the value of the input parameter of growth and 
density-dependence to the output for the density dependent model. Each plot shows 
the Pearson correlation coefficient and significance level: '.' p < 0.1 '*' p < 0.05 
'**' p < 0.01 '***' p < 0.001 
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Figure 11: Scatterplots relating the value of each input parameter of survival and 
fecundity to the output for the density dependent model. Each plot shows the 
Pearson correlation coefficient and significance level: '.' p < 0.1 '*' p < 0.05 '**' p 
< 0.01 '***' p < 0.001 
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Figure 12: Partial Rank Correlation Coefficients for the density independent (a) 
and density dependent (b) models. The bars are confidence intervals, generated by 
bootstraping 1000 times 
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Figure 13: eFAST analysis for the density independent (a) and density dependent 
(b) models. The bars represent the first and total order estimates for the sensitivity 
of each parameter in the model output. 
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variation caused by interaction between that parameter and others (total-order 
sensitivity). The interaction term for each parameter is the difference between its 
first and total order sensitivities. The eFAST analysis is substantially more intensive 
in terms of computer time than the previously mentioned. The analyses presented 
here required 7392 and 8448 runs of the simulations, respectively. Table [2] shows 
the SBMA measure of concordance between different sizes. Note that N s reported 
should be multiplied by the number of parameters in each model to obtain the total 
number of simulations executed. Notice that, while the main effect Di converge for 
both model, the total order D t indices are still variable with large sample sizes. 





Size (N s ) 


Indep. Di 


Indep. D t 


Dep. Di 


Dep. D t 


1 


66-132 


0.78 


0.85 


-0.14 


0.52 


2 


132-264 


0.83 


0.78 


-0.00 


0.75 


3 


264-528 


0.97 


0.93 


0.95 


0.66 



Table 2: Comparison of eFAST analyses by sample size for both models 



This analysis reveals that the output of the density independent model is mostly 
explained by first-order relations (which explain 62 % of the output variation), with 
F-j and the growth terms being the most important. The importance of the linear 
terms shouldn't come as a surprise, as the matrix growth model is a linear model. 
The density-dependent model, which has 66 % of the output variation predicted 
by linear terms, exhibit more complex interactions between the terms, in particular 
interactions involving s 7 , g m and a, the survival for adult plants, and the terms 
associated with competition between seedlings, respectively. 

4.3 Conclusions 

The results from the uncertainty and sensitivity analyses presented show some of the 
advantages from the methodology described in this work that are unavailable to the 
usual framework used in ecological studies. First, we have been able to quantify the 
uncertainty in the assymptotic growth rate (related to A) and the stable population 
size related to the uncertainty in the model inputs. Also, we provide a common 
framework to investigate side-by-side the linear and non-linear matrix models, from 
which we were able to point out the similarities and discrepancies between the 
models. Analyses based on matrix elasticities are also unable to investigate the role 
played by parameters not directly present on the matrix, as the size of the adult 
trees canopy z in our case. Finally, our approach permits the identification and 
quantifiation of relative importance of non-linearities and interactions between the 
input parameters in determining the model's outcome, and allows us to incorporate 
our previous knowledge about the system in specifying the range and distribution 
of each input parameter. 
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5 Previous use of Latin Hypercube in ecology 



The importance of sensitivity and uncertainty analyses in the development and 
use of ecological models is widely recognized. Searching for the terms "Sensitivity 
Analysis, Uncertainty Analysis or Parameter Space Exploration" in the Web of 
Knowledge reports 1199 papers in eight major journals (see fig. 14 and legend for 
details) since 1971, with 31650 total citations. However, most of these papers rely on 
full and individual parameter space exploration, which, as discussed in section [2j are 
not optimal. When restricting those results with the keywords "Latin Hypercube, 
MCMC, Markov, Monte Carlo", just 120 works show up in the results. Of those, 
only 13 (about 1%) use Latin Hypercube Sampling |31EII251I331H51I2SIESIIZII3DIE21 
[28| [T3l [21] . There are also relevant examples of LHS use in other journals [H |9| |4"2]. 

Also, many of these papers did not explicitly take into account the correlations 
between parameters. Those who did used mostly Iman and Conover's method |17j . 

These works have used Latin Hypercubes typically from 10 to 30 dimensions, 
but ranging from 6 to 143 [3], and the number of simulations ranged from 19 [28] 
to 2000 [43]. Also, these works are from varied areas within ecology: applied plant 
ecology [6j|43], species richness [13], epidemiology [36] and food chain analysis [7], 
stressing that the method is useful on varied problems. 
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Figure 14: Top: Number of papers per year since 1990 containing the topics "Sensi- 
tivity Analysis, Uncertainty Analysis or Parameter Space Exploration" in the jour- 
nals American Naturalist, Ecology, Journal of Ecology, Oikos, Oecologia, Ecological 
Modelling, Ecology Letters, Journal of Theoretical Biology as reported by Thomson 
Reuter's Web of Knowledge. Bot: Restriction of the search above to the keywords 
"Latin Hypercube". Search conducted 18.06.2012, 14h GMT. 
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